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Abstract 

To study and forecast the solar activity data a quite perspective method of singular 
spectrum analysis (SSA) is proposed. As known, data of the solar activity are usually 
presented via the Wolf numbers associated with the effective amount of the sunspots. 
The advantages and disadvantages of SSA are described by its application to the series 
of the Wolf numbers. It is shown that the SSA method provides a sufficiently high 
reliability in the description of the 11-year solar cycle. Moreover, this method is appro- 
priate for revealing more long cycles and forecasting the further solar activity during 
one and a half of 11-year cycle. 



1 Introduction 

It has long been observed that solar activity 
depends on a number of spots visible on its 
disk. During about 11 years which is called 
a solar cycle, this number varies over a wide 
range. Accompanying this process changes in 
the structure of magnetic fields of the Sun af- 
fect the Earth climate and have a probable 
connection with the natural catastrophes. In 
addition, intensity of the solar radiation (the 
frequency of the Sun flares, etc.) apparently 
exerts an influence on all of human life 
including socially-historical activity. Thus, 
in view of the quite large significance of the 
magnetic activity of the Sun its prediction is 
a subject of much current interest. 

At the present time, to describe the dy- 



namics of the solar activity many approaches 
are used. Among them the Wolf number 
associated with the effective amount of the 
sunspots is much more convenience method. 
The dynamics of the Wolf numbers has more 
or less quasiperiodic nature but, owing to the 
facts that models of this process do not take 
into account many essential factors of the so- 
lar magnetic activity, its prediction is diffi- 
cult. It is known that during last 250 years 
the cycle period of the sunspots changes its 
value not more than 20%. At the same time, 
the amplitude (i.e. the averaged number of 
sunspots) has varied over an order or even 
more. Detailed models of the solar processes 
do not describe such variations. 

In the last few years sufficiently many 
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methods devoted to the prediction and recon- 
struction (to the past) the dynamics of the 
Wolf numbers have been proposed (see, e.g., 
|1|, 0, 0, f|, U and references cited therein) . It 
should be noted however, that these methods 
have certain limitations and demerits. That 
is the reason why prediction of the sunspot 
dynamics based on the data of observations 
only (i.e. without a modelling of the process) 
is a quite perspective approach. In this way 
time-series analysis (see || [7|, |], || [10|, [II 



can give an essential assistance. But in this 
case there are many problems related to the 
fact that the Wolf series is not strictly de- 
terministic system, it does not have clearly 
defined dimension |E| 13|, and its length is 
not so large. 

In the present paper, for investigations 
of the Wolf time series a sufficiently efficient 
method of singular spectrum analysis is pro- 
posed. It is shown that this method provides 
a quite high reliability in the description of 
amplitudes of the 11-year solar cycle. More- 
over, it is appropriate for revealing more long 
cycles and forecasting the further solar activ- 
ity during one and a half of 11-year cycle, i.e. 
the nearest 16 years. 

2 Singular spectrum 
analysis 

The method of singular spectrum analysis 
(SSA) 0, |15|, U, 0, |U, which is used in 
the present article allows us the following. 

• Recognise certain components in the 



time series which have been obtained 
from the observable at regular intervals; 

• Find periodicities that are not known 
in advance; 

• On the basis of the chosen components 
to smooth out the initial data; 

• Extract components with the known 
period; 

• Predict the further evolution of the ob- 
served dependence. 

The method of SSA is a reasonably new 
one, but now it is quite clear that it is suf- 
ficiently competitive with numerous smooth- 



ing methods |L9], ^0], [ZI] . Moreover, in certain 
cases forecasting the system evolution on its 
basis gives much more reliable results in com- 
parison with the other known algorithms (see 
f7|, 118, [22], [24], |25| and references therein) . 

The SSA method is based on the pas- 
sage from the investigation of an initial lin- 
ear series (xi)f =1 to the analysis of a many- 
dimensional series consisted of components 
of some length which contain except for the 
value Xi, certain quantities x^j, j = 1, . . . , r, 
at the previous instants of time. 

Let us describe the central steps of the 
application of SSA to a series {xi)f =l . 

(1) On the first step, a one-dimensional se- 
ries is transformed into a many-dimensional 
one. For such a transformation it is nec- 
essary to take a certain number of delays 
t < W/2 + 1], where [•] is an integer part, 
and represent the initial r values of the se- 
quence as the first column of some matrix X. 



2 



For the second column, the values of the se- 
quence from x 2 to x T+ i are chosen. Thus, 
the last t elements of the sequence x n , . . . , xn 
correspond to the last column with the num- 
ber n = N—t+1. Therefore the transformed 
series has the following matrix form: 
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trix of the conversion to the principal com- 
ponents, Y = V T X = (Y u Y 2 ,..., Y T ), of the 
initial series. Here Y^, i — 1, 2, . . . , r, are rows 
of the length n. Therewith, the eigenvalues 
Ai, A2, . . . , A T can be considered certain 
contribution of the principal components to 
a general information content of the series 
{xi)f =1 . Then, by means of these principal 
components it is possible to reconstruct the 
initial matrix X: 



J 

It is obvious that for the constructed matrix 
the expression \\xij\\ = Xi + j-\ holds. In gen- 
eral, the matrix X is a rectangular one. But, 
in a limit case, i.e. at r = N/2 and an even 
N, it degenerates into a square matrix. 

(2) After this transformation, for the ma- 
trix X the corresponding covariance matrix 

1 



x = (v\v 2 ,...,v T 



Y 2 T 
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should be obtained. 

(3) Now it is necessary to find eigenvalues 

and eigenvectors of the matrix C. For this, 

the matrix C should be factored as follows: 

/ Ai ... \ 

A 2 ... 

C = VAV T , where A = 

V ... K J 

is the diagonal matrix of eigenvalues, V 



In turn, by the matrix X one can recon- 
struct the time series (xj)^. It should be 
noted that for the reconstruction not all prin- 
cipal components Y 1 ,Y 2 ,...,Y r are usually 
applied. Only a part of them can be involved. 
This depends on the goal which we pursue 
and the informative content of the used com- 
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the orthogonal matrix of eigenvectors of the 
matrix C. It is clear that A = VCV T , 
Y7i=\ Aj = t and det C = nj=i A«. 

(4) For the next step, the matrix V of 
eigenvectors should be presented as a ma- 



ponents (see [0, [15], [L6|, [T7|). This means 
that each vector-row Yi can be considered as 
a result of some projection of a r-dimensional 
totality on a direction corresponding to the 
eigenvector V\ Thus, the series is presented 
via a set of r components Yi. Therewith, the 
weight of the component Y^ in the initial se- 
quence (xi)f =1 can be defined by the corre- 
sponding eigenvalue A» which is, in turn, the 
eigenvalue of the eigenvector V 1 . 

Each i-th eigenvector includes r compo- 
( v\\ 



nents, V 1 



Let us construct a de- 



V vl J 



pendence of the component values v l k , k 
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1, 2, . . . , r, as a function of their number: 
v l = v l (k). Then using the orthogonality 
property of eigenvectors, the further analy- 
sis of the sequence (xi)f =1 can be performed 
by means of diagrams plotted by the anal- 
ogy with the Lissajous figure. Namely, along 
the axes the components v l k , v{ are plotted in 
pairs. If the constructed diagrams are simi- 
lar to a circle, then the functions v l = v l (k), 
i)i = v J (k) are approximated by certain peri- 
odic functions with almost coincided ampli- 
tudes and a phase lag about a quarter of the 
period. 

Thus, for some pairs of the eigenvectors 
V 1 , V 3 one can find a value which has a mean- 
ing of the period. Therefore, such a graphi- 
cal representation provides a certain estimate 
of the component frequencies inherent in the 
initial time series {xi)f =1 . 

For the given parameter r the number of 
all possible pairs of principal components is 
~ t 2 . It is obvious that even for a suffi- 
ciently small t the analysis of all such pairs 
is a quite laborious problem. Moreover, at a 
large values of r only a small part of diagrams 
has a helical form. Thus, before a graphical 
analysis it would be reasonable to restrict our 
search. This can be done if we arrange V 1 and 
Yi in order of decreasing their eigenvalues and 
consider only such pairs of eigenvectors which 
have close enough eigenvalues A. In the dia- 
gram A = at a quite large A these pairs 
form a decreasing (with the growth of i) step 
function. 

(5) Suppose now that for the further re- 
construction we have only r leading compo- 



nents. Thus, for the reconstruction of the ini- 
tial matrix X one should use r leading eigen- 
vectors V 1 . In this case 



x = (v\v 2 ,...,v r 



( n \ 

\ Y r J 



i=l 



where X is a reconstructed matrix with n 
columns and r rows. Then the initial time 
series obtained from this matrix is defined as 
follows: 



f 1 s 
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1 < S < T , 

t < s < n , 
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7l< S < N 



Described way of the reconstruction is said to 
be a SSA-smoothing of the initial time series 



1 X i 



1 by the leading r components. 



(6) In the next stage of the SSA applica- 
tion, a prediction procedure of the initial time 
sequence can be considered (see 
This means that the series 
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x i)i=L P which 



is the extension of the known data 



, x 



i)i=li 



is constructed. In turn, extrapolation to p 
points forward is reduced to the application 
of p times of the prediction procedure to the 
one point. The basic idea of the computation 
of the point xjy + i is the following. 

Consider the sequence xi, x 2 , ■ ■ ■ , xn and 
construct a sample in the form of matrix X. 
As a basis of the surface containing this sam- 
ple one can take the chosen before vectors 
V 1 , V 2 , . . . , V r of the matrix C. Let us write 
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the parametric equation of this surface as 

T 

S(P) = PiV' \ where a set of the param- 
i=i 

eters Pi corresponds to the value S(P) which 
is a column with r elements. In this case the 
set of parameter values P k = (pi,p^, ■ ■ ■ ,p^) 
corresponds to fc-th, k = 1,2, ... ,n, column 
of the matrix X. Therefore, X 1 = ^(P 1 ), 
X 2 = S (P 2 ), ...,X n = S (P n ). Now, to pre- 
dict the value xjy + i it is necessary to find the 
(n+l)-th column X n+l which, in turn, fits the 
parameters P n+1 = (pi +1 ,p% +1 , . . . 
Using the data {xi)f =1 these parameter values 
can be obtained from the expression S(P) = 

J2 PiV 1 . Thus, the predicting column is writ- 

i=i 

ten as follows: X N+1 = S(P n+1 ). 

Let us introduce the following designa- 
tions: 







( v\ 


v 2 . 
v 2 . 


■ Vl \ 

. vl 










■ <-! ) 




p = 


( Pl +1 \ 

P n 2 +1 


; Q = 


( XN-t+2 
XN-t+3 


\ 




\ P n r +1 ) 




V X N 


/ 




v T = 


Uy T ,... 


•<)■ 





The set of the values (pl +l ,p r 2 l+l , . . . ,p™ +1 ) 
is easily found solution of the system 

V*P = Q with respect to P. Thus, the final 
expression of the forecasting value has the fol- 
lowing form: 

KK T Q 

In the simplest case, to predict the next 
values it is necessary to change the matrix Q 



in a corresponding way and multiply it by the 
value V T V? / (l - V T V T ). In addition, how- 
ever for each predicted point one can com- 
pletely repeat the SSA algorithm. Then the 
matrixes V T and V* will be changed. 

(7) At the final stage of the SSA ap- 
plication one should dwell on the choice of 
the main parameter — the number of delays 
r which are used for the many-dimensional 
sample X. As in the case of selection of the 
principal components this value essentially 
depends on the investigated problem. 

Consider the smoothing procedure of a 
time series by the SSA method. In this case, 
as noted above, selection of a principal com- 
ponent is the filtration of the time series with 
the transition filtering function in the form of 
the eigenvector of this principal component. 
If the delay value r is the greater, the greater 
the number of parallel filters, and a band- 
width of each of them is more narrow. Thus, 
for a large enough r we have a sufficiently 
efficient smoothing of the time series. 

If it is necessary to define unknown (hid- 
den) periods in the observed sequence then 
one should take the value of r as large as 
possible. Next, after omitting close to zero 
eigenvalues, the delay value should be short- 
ened. 

Suppose now that we should select only 
the one known periodic component. In this 
case it is necessary to choose the delay r 
which is equal to the required period. 

Finally, let us consider the problem of 
some extension of an observed sequence to 
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a given value, i.e. the problem of forecast- 
ing the evolution of the process under inves- 
tigation. Then one should take a maximum 
allowable value of the delay r and thereafter 
select the number r. 

3 Forecasting the solar 
activity by SSA 

To estimate of the practicality of the SSA 
method it is necessary to use the time se- 
ries of a natural origin. In the present ar- 
ticle, as such a time series the sequence of 
the Wolf numbers characterised in a certain 
way the solar activity was chosen. Taken 
alone, the Wolf numbers defined by the visi- 
ble sunspots cannot give a quantitative infor- 
mation related to the solar activity. However 
there is a large enough correlation between 
the Wolf numbers and the F10.3 emission. 
That is the reason why investigations of the 
relative variations in this sequence can give a 
specific information concerning the solar ac- 
tivity. 

For the first time, in 1848 a Swiss as- 
tronomer R.Wolf proposed that a measure 
of the solar activity can be characterised by 
the number of sunspots. To this end he rec- 
ommended to consider the union of the total 
number of spots visible on the Sun and ten- 
fold numbers of regions in which these spots 
are placed. This last summand should coor- 
dinate the results of measurements performed 
under several conditions. Thus, since the 
year 1849 the results of daily measurements 
have become available. Using more earlier 



observations and various sources, R.Wolf has 
reconstructed the data of the solar activity 
(with an admissible accuracy and negligible 
gaps) up to the year 1818. Now the averaged 
number of the sunspots is called the Wolf 
number. 

Later the average monthly values of the 
Wolf numbers up to the year 1749 (namely 
this series is used in the present paper) and 
their average annual values up to the year 
1700 have been reconstructed. In the last 
case however, the error can reach more than 
10 percent. The chosen data covers the wide 
time interval without gaps and a quite high 
time resolution. The investigated sequence is 
shown in Fig.l. The time with the interval 
of one month is plotted on the abscissa. The 
corresponding Wolf numbers are laid as ordi- 
nates. 

On the first step of the application of SSA 
it is necessary to take the maximum permis- 
sible value of parameter r. For our investiga- 
tions we used r = 500. This value allows us 
to make up periodicities up to period of 42 
years. The use of the larger values leads to 
essential computer problems. Moreover, in- 
creasing the parameter r (up to 600) does not 
yield any significant change in the results of 
the first expansion in principal components. 

Owing to a large value of r the depen- 
dence of the roots of eigenvalues of the co- 
variance matrix (ordered in decreasing) de- 
cays exponentially. In addition, the sum of 
the leading five eigenvalues exceeds 99% of 
their sum total. Coupled with a quite large 
numbers of initial points, this leads to the 
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fact that the first principal component yields 
a small enough smoothing of the initial series, 
and by means of the leading four-five compo- 
nents one can reconstruct this series. More- 
over, at small values of the parameter r, say 
at r = 5, the first component weakly changes 
its form. This is due to the fact that the SSA 
method is stable with respect to this param- 
eter. Therefore, application of a sufficiently 
large r is justified only for the prediction. 

To check possibilities of the prediction let 
us cut off the initial sequence of the aver- 
aged monthly values for a length of 216 points 
(18 years) and reconstruct it by the follow- 
ing way. Evaluate optimal parameter values 
of the reconstruction algorithm by means of 
an additional truncation and resolve the ob- 
tained series into r components. Therewith 
it is necessary to choose such a number of 
the leading components at which we get the 
best coincidence of the reconstructed values 
with the additionally truncated data. Then 
using the obtained parameters, we will find 
the initially shortened part of 216 points. 

By means of the direct exhaustion one 
can define that the best results we get at 
r = 150 (the number of the chosen compo- 
nents). Once again let us take the reduced 
sequence (for a length of 18 years) and apply 
the chosen r for the prediction. One can fur- 
ther improve the quality of the prediction if 
we will decompose the predicted interval (216 
points) into subintervals and after prediction 
of each part calculate the principal compo- 
nents. In the ideal case, after the prediction 
of each point it is necessary to realise this 



process. However, such an approach leads to 
an increase of the computer time. 

Fig.2. illustrates the reconstructed series 
at which recalculations have been realised ev- 
ery 72 points (i.e. 3 times). However, as it 
turns out, this result is almost identical to the 
prediction of the whole interval of 216 points, 
without decomposition into subintervals. 

Theoretically, one can analyse the series 
components with the aim of detection the ex- 
istence of the known and any other periods. 
However a number of these periods and thus, 
the resemblance of the components Yi of the 
initial series makes this analysis a quite cum- 
bersome procedure. In addition, selection of 
sufficiently short periodicities (about several 
months) is a very difficult computer task. For 
these reasons, the choosing a large enough 
subinterval of time is the more simply prob- 
lem. 

Let us consider now a series with the av- 
erage annual Wolf numbers. This series con- 
tains only 248 points. Thus, the maximum 
allowable value of r is 123. Choose this value 
as an initial one. The corresponding eigen- 
values is shown in Fig.3. The first number 
presents the principal component related to a 
trend. The following step-like data form the 
pairs of components with the numbers 2-3, 4- 
5, 6-7, 8-9 and 11-12. Beginning with 14-th 
number this dependence gives way to an ex- 
ponential tail. Eigenvectors for the pairs 2-3, 
4-5, 8-9 and 11-12 (see Fig.4) fit periodici- 
ties with 11-year period. The corresponding 
helical dependence for the components 2 and 
3 is shown in the left side of Fig. 5. 
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It should be noted that except for the ev- 
ident the eleven-year solar cycle it is possible 
to guess a supposed eighty- year Gleisberg cy- 
cle (see Fig. 5). Here we keep in mind the pair 
of the eigenvectors 6 and 7. The correspond- 
ing eigenvalues is not exactly equal to each 
other, i.e. the step-like data is oblique (see 
Fig.3) and the phase lag is not ir/2. That is 
the reason why the form of the function dia- 
gram is not helical. In spite of this fact and a 
small enough eigenvalues, it is quite possible 
to trace a periodicity. 

For the best selection of the eighty-year 
dynamics one can identify parameter r. By 
the numerical analysis we have found that 
r = 80 is the most suitable value. In this 
case, vectors 4 and 7 are associated with such 
a periodicity. The diagram for 4-th and 7- 
th components is shown in Fig. 6. A possible 
eighty- year solar cycle obtained by the recon- 
struction only via these two components is 
shown in Fig. 7. 

Consider now a much more interesting 
problem concerning the possibility of the pre- 
diction of the average annual Wolf numbers 
by the SSA method. Shorten this series from 
the right hand side for a length of 18 points 
(that means years) and reconstruct it. To 
make this let us truncate additionally an in- 
terval of 11 points and try to restore it in the 
best way. Choosing the suited numbers of 
eigenvectors we will use such a procedure in 
various parts of the series. 

For the series of 219 points (219 = 248 - 
18 — 11, where 248 is the number of all points, 
see above) the maximal possible r = 109. As 



follows from numerical analysis, prediction 
for 11 points strongly depends on the pick- 
ing of the components. At the same time, in 
a sufficiently wide range a qualitative guide 
of the prediction is not bad. However, for 
the quantitative prediction the given r is too 
large. For the lesser values of r it is quite pos- 
sible to find the necessary number of compo- 
nents. For example, for r = 33 and choosing 
the leading 11-th components we get fairly 
good results of the prediction. 

Thus, let us use these parameters for the 
prediction of the removed 18 points. In con- 
trast to the case of the average monthly data 
and r = 500 now, at r = 33 we have a pos- 
sibility to recalculate eigenvectors with con- 
sideration of the last predicted point. This 
recalculations can be made at the stage of 
the selection of vectors as well as the predic- 
tion. The result of the prediction with the 
correction at each step is shown in Fig. 8. 

In addition, we have studied a sequence of 
natural logarithms obtained from the initial 
series. Taking the logarithm is often used at 
the data processing (for example, in the case 
of the correlation analysis) and allows to get 
more interesting and complete results. How- 
ever, eigenvalues and eigenvectors after such 
an operation have not principally changed. 
Thus, in this case taking the logarithm is not 
necessary because the basic information can 
be obtained from the analysis of the initial 
series. 

In the closing stage let us consider appli- 
cation of SSA for a real prediction of the solar 
activity. For this purpose, the average annual 
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sequence of the Wolf numbers from the year 
1748 to the year 1996 has been chosen. The 
end of this series corresponds to a minimum 
of the solar activity. Therefore it is interest- 
ing to describe the further activity of the Sun 
and predict its two next maxima. To realise 
this idea it is necessary to extend the aver- 
age annual sequence for a length of 18 year 
points. 

We resolved the series into 33 compo- 
nents and choose for the prediction the lead- 
ing eleven of them. The result of the predic- 
tion is shown in Fig.9. As follows from this 
figure, in the nearest future, in comparison 
with the two previous maxima, the Sun will 
be in a relatively quiescent state. In addition, 
the level of the forthcoming maximum in the 
year 2011 will not be so high. On the ba- 
sis of our numerical investigations concerning 
the application of SSA, the prediction during 
more than one and a half of 11-th year solar 
cycle is not justified. Nevertheless, we believe 
that the obtained values of the Wolf numbers 
shown in Fig.9 are a quite reliable. 

For comparison, a series of the average 
monthly Wolf numbers has been analysed. 
However, its investigation did not yield an es- 
sential modifications to the obtained results. 

4 Conclusion 

It seems to be true that in the nearest fu- 
ture investigations of time series by means of 
SSA will occupy a more important and de- 
serving place among different ways related to 
the processing and forecasting many experi- 



mentally obtained sequences. Resolving the 
initial series into components the analytical 
form of which is not fixed, this method allows 
us on a quite high level to reveal periodic sub- 
sequences and forecast the dynamics of this 
series. Therewith, restrictions to the num- 
bers of points and characteristic periods in 
the investigated data rule, is essentially 
less than at applications of the other methods 
(for example, at correlation or Fourier analy- 
ses). 

In the present article a possibility of the 
SSA application to the sequence of the Wolf 
numbers characterising the solar activity, is 
considered. In spite of a relatively small 
length of this sequence, SSA allows to reveal 
the components corresponding the known so- 
lar cycles and, on the basis of only certain 
constituents, admits reconstruction of the 
data. Moreover, it is found that by means 
of SSA it is possible to extend short time se- 
ries with an acceptable accuracy. 

Like any other method, SSA is not de- 
void of drawbacks. First, there is a cer- 
tain difficulty in the problem of revealing 
unknown frequencies in the investigated se- 
quence. These frequencies can be more easily 
obtained by the Fourier analysis. Secondly, 
SSA does not contain clear rules concerning 
the choice of components, especially in the 
case of forecasting. Finally, SSA does not 
give the correct prediction of the cycle period. 
Therefore in the predicted sequence a system- 
atic phase lag is accumulated. Nevertheless, 
as follows from the performed analysis it is a 
useful supplement to the existing methods of 
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the experimental data processing. 

In the forthcoming investigations devoted 
to the analysis of the solar activity, we will 
consider a possibility of a certain correction of 
the cycle phase by means of additional meth- 
ods and via special feature of the sunspot dy- 
namics. Besides, to resolve this problem one 
can use an empirical dependence of the am- 



plitude and the phase of the solar cycle |PS|| . 
Moreover, it seems to be reasonable that it is 
quite possible to improve control parameters 
which are used in SSA. 

On the whole we can say that the de- 
scribed method of singular spectrum analysis 
is a sufficiently advantageous and promising 
way for the prediction of dynamics of the so- 
lar activity. 
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Observations — □ — 




Figure 2: An example of the prediction of the solar activity for a length of 216 points (18 
years) by the average monthly Wolf numbers. The vertical line corresponds to the boundary 
value of the removed points. At the expansion, 500 components have been applied; in the 
reconstruction procedure 150 components have been involved. Numerical analysis has been 
performed by three stage: after prediction of the next 72 points recalculations have been 
made. 
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Figure 3: The leading 50 eigenvalues of the covariance matrix obtained at the expansion of 
the average annual Wolf numbers into 123 components. 
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Figure 4: Some of components corresponding to eigenvalues shown in Fig.3. Their percent- 
age in the initial series is denoted by the parentheses. 
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Figure 5: Dependencies of the 2-nd and the 3-rd components (on the left-hand side) and 
the 6-th and the 7-th components (on the right-hand side) shown in Fig. 4. The left and the 
right parts of the Figure correspond to the 1-st and the 3-rd step-like data of eigenvalues, 
respectively. 
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Figure 6: Dependencies of the 2-nd and the 3-rd components (on the left-hand side) and 
the 4-th and the 7-th components (on the right-hand side) at r = 80. 
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Figure 7: Reconstruction of the 4-th and the 7-th components corresponding 80-year solar 
cycle. 
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Figure 8: Prediction of the solar activity the each-step correction. 11 in 33 components are 
chosen. The vertical line marks the boundary value of the removed 18 points. 
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